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ABSTRACT 

Characterization of migration in gravitationally unstable disks is necessary 
to understand the fate of protoplanets formed by disk instabihty. As part of 
a larger study, we are using a 3D radiative hydrodynamics code to investigate 
how an embedded gas giant planet interacts with a gas disk that undergoes 
gravitational instabilities (GIs). This Letter presents results from simulations 
with a Jupiter-mass planet placed in orbit at 25 AU within a 0.14 Mq disk. The 
disk spans 5 to 40 AU around a 1 Mq star and is initially marginally unstable. In 
one simulation, the planet is inserted prior to the eruption of GIs; in another, it is 
inserted only after the disk has settled into a quasi-steady Gl-active state, where 
heating by GIs roughly balances radiative coohng. When the planet is present 
from the beginning, its own wake stimulates growth of a particular global mode 
with which it strongly interacts, and the planet plunges inward six AU in about 
10^ years. In both cases with embedded planets, there are times when the planet's 
radial motion is slow and varies in direction. At other times, when the planet 
appears to be interacting with strong spiral modes, migration both inward and 
outward can be relatively rapid, covering several AUs over hundreds of years. 
Migration in both cases appears to stall near the inner Lindblad resonance of a 
dominant low-order mode. Planet orbit eccentricities fluctuate rapidly between 
about 0.02 to 0.1 throughout the Gl-active phases of the simulations. 

Subject headings: hydrodynamics — instabilities — planet-disk interactions — planets 
and satellites: formation — protoplanetary disks 
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Introduction 



Migration of gas giant planets due to interactions with a circumstellar gas disk can 
play a major r ole in defining the arch itecture of planetary systems. Work on migration 



(see review by 



Papaloizou et al 



20071 ) has included gravitational interaction of planets 
with both laminar and turbulent disks. However, radiative transport, detailed equations 
of state (EOS), and the self-gravity of the gas disk have usually been ignored; the effects 



of a non-isothermal EOS 



2006 



lave only re cently been included (e.g.. 



Paardekooper et al. 



2010 



Paardekooper fc Mellema 



201 ll ). Emerging studies, such as 



Baruteau et a. 



(2011) and the one repo r ted h ere, are beginning to address some of these issues. 



Boss 



(120051) and 



Mayer et al. 



( 120041 ) examined radial migration of planet-mass fragments in 
gravitationally unstable disks, but their disks were violently disrupted by fragmentation 
under conditions (radii < 40 AU, disk-to-star mass ratios M^/ M.. ~ 0.1, and stellar 



mass M.. ~ 1 Mp-.) where fragmentation may not actually occur f Rafikov 



Boley et al. 



2006 



2007a 



Boley Sz Durisen 



2008 



Forgan et al 



2009 



2005 



Cai et al. 



2007 



2010|). More 



recently, fragmentation into clumps with gas giant or brown dwarf masses has been 
documented in numerical simulations of disks that are relativ ely massive (Md/M, 



a few tenths) and spat i ally extended (outer radii > 50 



Stamatellos et al. 



2007 



Stamatellos fc Whitworth 



2009 



AU) 



Boley 



( Krumholz et al. 



2007 



20091: 



wher e fragmentation is expected fr om semi-analytic arguments (e.g. 



2009 



Dodson-Robinson et al. 



Bolev et al. 



Clarke 



2009 



2010h. 



Rafikov 



20091 ) . The fate of the clumps then depends in part on the ir 



radia l migration, which is a chaotic and messy affair in a fragmenting d isk (e.g. 



2009 



Boley et al 



2010 



Vorobyov &: Basu 



2010a 



Bolev &: Durisen 



of gravita 


;iona 


Zhu et al. 


2010) 



instabilities (GIs) may be episodic (e.g.. 



2010) 



Boley 



Vorobyov fc Basu 



'he occurrence 



20061 . 



2010b 



later find themselves in a disk that erupts again into GI activity. As the star /disk system 
evolves, such a protoplanet may end up in a region of a Gl-active disk where fragmentation 
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does not occur. 



To improve our understanding of how planets migrate in Gl-active disks that 
are not fragmenting, we have begun a systematic study, using numerical 3D radiative 
hydrodynamics, where we investigate the effects on both the disk and the planet of 



inserting a planet-mass object into disks susceptible to GIs. Using 



earlier research by our g r oup (iPickett et al 



Boley et al. 



2006 



2007a 



Michael fc Durisen 



2003 



Mejia et al. 



t echniques deve. 



2005 



Cai et al 



2006 



oped i n 



2008 



2010l ). we can identify the dominant spiral 



waves in a simulation and analyze how the waves interact with the planet's motion. Our 
goal is to determine both the effect of giant planets on GIs and the effect of GIs on planet 



migratio n. Because GIs ar e sensitive to radiative physics, we use a wel 



scheme (IBoley et al. 



2007al ) and realistic opacities (ID'Alessio et al.ll200ll ). 



-tested radiative 



Section 2 below presents our numerical methods and initial conditions. We describe 
the simulation results in S3, and discuss them in S4. 



2. Computational Methodology 



2.1. 3-D Radiative Hydrodynamics 



The CHYME RA (Computation al HYdrodynamics with MultiplE Radiation 



Algorithms) code ( IBoley et al. 



2007al ) is a second-order, explicit. En 



erian scheme on a 3D 



2007b[ ) 



cylindrical grid. The code uses a realistic equation of state for H2 ( IBoley et al. 
and integrates an energy equation that includes PdV work, net heating or cooling due to 
radiative flux divergence, and heating by artificial bulk viscosity. Calculations are done 
on a uniform cylindrical grid with reflection symmetry about the disk midplane and a 
grid size {w,(f),z) = (512,512,64). The z-axis is the rotation axis of the disk. The large 



number of azimuthal zones is necessary to resolve the planet's Hill sphere and the planet's 
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wake. These simu^ 



Boley et al. 



ations utilize the radiative coohng scheme developed and tested in 
( l2007al ). where the optically thick monochromatic flux in the w and 0-directions 
is computed by flux-limited radiative diffusion and where the radiative transport of energy 
in the 2;-direction is solved using a one-ray discrete ordinate method in both optically thin 
and thick regions. Although the central star remains flxed at the grid center, we account for 



acceleration of the re 



m 



Michae. 



Hut et al 



'e rence frame by the planet and by the disk via indirect potentials, as 



fc DurisenI (120101 ). The planet integration is done with a Verlet integrator (e.g.. 



19951 ). and the indirect potential terms are treated as in 



Nelson et al. 



mm . The 



Rosseland mean and Planck mean opacities an d molecular weights in our simulations are 



the same as those in 



Boley et al. 



(2006 



2007al ). except that we correct the 



D'Alessio et al. 



(I2OOII ) gas mean molecular weight to n = 2.33. 



2.2. Initial Model and Simulation Conditions 



Pickett et al. 



(120031 ). orbits a 



The model disk, based on an equlibrium model from 
1 Mq star and has a mass Md = O.14M0, inner and outer radii at 5 and 40 AU, and 
an initial surface density distribution S ~ -07"^/^. The time unit of one ORP (= outer 
rotation period) is deflned as the rotation period of the initial disk at m ^ 32 AU, or 
about 180 yr. The disk is initially located between radial grid zones 30 and 240 and is 
close to isentropic, which results in a Toomre-Q distribution with a marginally unstable 



see 



Durisen et al. 



20071 ) minimum value of 1.38 at radial grid zone 161 (26.7 AU). The 



computational grid extends radially to 512 zones to accommodate expansion of the outer 
disk once GIs become nonlinear. An outflow boundary condition is enforced at the upper 
vertical grid boundary, the outer radial grid boundary, and an inner radial boundary at 2 
AU. To seed nonaxisymmetry, the density distribution is given an initial 0.01 % random 
cell-to-cell perturbation. 
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This Letter presents three simulations. The first, which we call the Gducial run, is 
simulated from time t = to 21 ORP without a planet. In the second simulation, which we 
call the t = planet run, a one Jupiter-mass (Mj) planet with a roughly circular midplane 
orbit is included in the initial t = ORP equilibrium disk at 25 AU. In the third simulation, 
which we call the t = 10 planet run, a 1 Mj planet with a similar orbit is inserted into 
the fiducial run at t = 10.5 ORP. These three simulations are part of a larger suite, to be 
described elsewhere, in which the planet mass is varied. We estimate semi-major axes a and 
eccentricities e from the planet's w{t) by using the maximum and minimum radii, Wmax 
and ^min, over each complete orbit of the planet, where one orbit is = 27r. Specifically, 
a = (ti7max + Wmin)/2 and e = (winax " '^mhi) / {'^ma.x + ^^min)- Thcsc parameters are intended 
to highlight the strong orbit-to-orbit variations in the motion of the planet. A more refined 
estimate of the orbital elements is unnecessary because the potential of the disk, even if 
were axisymmetric, will cause departures from Keplerian dynamics. 



Results 



The Fiducial Run. As in iMejia et al.l ( 120051 ). the fiducial run exhibits four distinct 
phases: initial axisymmetric cooling, the onset of GIs in a well-defined burst of a few 
discrete low-order modes, a transition to more complex nonaxisymmetric structure, and the 
asymptotic establishment of quasi-steady GI activity with an overall balance of heating 



and cooling. There is no fragm entation because this disk has a long cooling time (IGammie 



2001; 


Bolev et al. 


2006. 


2007a) 



The t = Planet Run. Figure [T] shows that inclusion of even this modest mass planet 
(0.7% of Mrf) has a dramatic effect on the burst phase. Without a planet, when all modes 
grow from imposed noise, it takes 4 ORPs for coherent spiral modes to organize and grow 
to nonlinear amplitude. The growth is centered near the Q-minimum at 26 AU. Modes of 
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cos(m(/)) symmetry with m = 4 and 5 dominate. An m = 3 (three-armed) spiral also grows, 
but somewhat more slowly, and lags the other modes in time. In the t = planet run, the 
planet develops an organized wake within about an ORP in which m = 3 is a dominant 
component (~ 5 % global density perturbation). This nonlinear seeding causes m = 3 to 
dominate the GI burst, which also occurs ~ 2 ORP earlier. Because of our initial placement 
of the planet inside but close to the Q-minimum, the corotation radius (CR) of the m = 3 
mode is fairly close to that of the planet's orbit radius. When the triggered m = 3 mode 
becomes strongly nonlinear at about 3 ORP, planet migration is significantly affected. 

Figure [2] shows the evolution of the planet's radial position. From to 2.5 ORP, the 
planet is torqued primarily by its own wake and migrates inward. Beginning at about 3 
ORP, the planet interacts with the now nonlinear m = 3 GI mode. At first, the planet gains 
angular momentum and moves outward, but, from 4 to 8 ORP, a time interval of about 720 
yr, the plant experiences a negative torque and plunges from 23 to 17 AU. After t = 8 ORP, 
the main burst is over, and the disk transitions into its asymptotic state, where modes 
of many m- values become comparably strong (see Fig. [1]). The planet's radial migration 
apparently stalls at about 16 to 17 AU. From an analysis of periodicities present in the gas 
disk between 8 and 12 ORP, the planet lies a few AU inside the inner Lindblad resonance 
of a strong m = 2 mode with CR at 29 AU in this early part of the asymptotic phase. 

The t = 10 Planet Run. The motion of the planet in this case is more difficult to 
interpret. We see intervals of fairly rapid inward or outward migration over several orbits, 
as well as times when radial migration appears to stall. Between t = 15 and 19 ORP, the 
pattern of outward migration for 2 ORP followed by an inward plunge over the next 2 ORP 
resembles the behavior in the t = planet run during the burst, but, in this case, there 
is no distinct transition between phases of GI activity. Animations of the evolution of the 



midplane density (available at 'http: / /hdl. handle.net /2022 / 13304 ) show there is a complex 
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interaction between the planet, its spiral wake, and the global spiral arms caused by the 
instabilities. Analysis of images and periodicities suggests that, between 15 and 19 ORP, 
the planet may be in a mean motion resonance with m = 2 and 3 patterns, both of which 
have CR at 27 AU. Over this time interval, the average orbit period of the planet is ~ 0.50 
ORP, while both the m = 2 and 3 spirals have a pattern period of ~ 0.74 ORP, suggesting 
a 3:2 resonance. Both of these modes are present, at nearly the same periodicities, in the 
fiducial run. So, in contrast with the t = planet run, this planet's rapid phases of radial 
migration are caused by interaction with strong GI modes that exist independent of the 
planet. 

Eccentricity. Figure [3] shows the evolution of e in the planet simulations. In both 
cases, the planets were inserted with an approximate circular velocity computed by adding 
the interior disk mass to the stellar point mass. The presence of GI activity in the t = 10 
planet run immediately increases e to ~ 0.08. In the t = planet run, the modest initial e 
decreases during the 1.5 ORP when it is naigrati ng only due to its own wake, as expected 



( IWard fc Hahn 



1998 



2003 



Goldreich et al. 



2004J ). However, once there are strong nonlinear 



interactions with GI modes, e jumps upward. In both runs, it appears that interaction 
with well-established GI activity leads to eccentricities between 0.02 and 0.1 that vary in 
a chaotic way between these extremes on orbit period time scales. The magnitude of e is 
roughly consistent with the ratio of the sound speed to orbital speed of the disk and the 



modest values for the Mach numbers of the GI spirals (iBoley fc Durisen 



2006 



20081). In 



other words, the planet's orbit tends to have about as much non axisymmetry as refiected 



in the pitch angles of the Gl-induced spirals (iCossins et al. 



20091). 



Migration Relative to the Disk. The disk surface density distribution does not vary 
much in the asymptotic phase, so radial motion of the planet in the t = 10 planet run a 



represents radial migration relative to the gas. During the burst (e.g.. 



Boley et al. 



20061), 
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the surface density of the disk changes dramatically. To verify that the planet is really 
migrating relative to the background gas disk in the t = run, we have compared the 
evolution of the planet's radial position Ci7p(t) with the mass inside that radius, expressed 
in terms of rricyi, the fraction of the disk mass interior to the planet. Between t = 3.52 and 
3.89 ORP, when tUp increases modestly, rricyi increases from about 0.50 to 0.53; between 
3.89 and 8.47 ORP, when Wp decreases dramatically, rricyi also decreases dramatically from 
0.53 to 0.30. So, the planet is migrating significantly with respect to the disk mass even 
while the inner part of the gas disk is moving inward during the burst. As a measure of 
the gas motion, rricyi at 23 AU increases from 0.50 to 0.60 from t = 3.89 to 8.47 ORP, 
corresponding to an average gas disk inflow rate of ~ 1O~^M0 yr~^. 



4. Discussion 



Even in non-self- gravitating disks, the presence of discrete radial str ucture is known 



(Li et al. 


2000. 


2001 





200ll ). iMeschiari fc Laughliru (120081 ) demonstrated 



that a gap, presumably opened by an embedded planet, can cause global modes to become 
unstable in self-gravitating disks that are otherwise stable against the development of GIs. 
Here we are dealing with a planet that is not quite massive enough to open a gap, but its 
gravitational interaction with the disk also has a strong effect on the growth and onset of 
GIs. 



Figure [2] suggests that the overall radial drift during the first 10 ORP of the t = 
planet run is not too different from what one would get by extrapolating the migration due 
solely to the planet's wake during the first 2 or 3 ORP, where the average radial migration 
is about 4x10^^ AU yr^^. Although there are no s imple formula e that are fully applicable 



to our radiatively cooled disk, using eq. (70) from 



Tanaka et al. 



(I2OO2I ) with parameters 



from the simulated disk, we estimate 6x10 AU yr ^ as an analytic expectation. This is a 
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reasonable estimate for compariso n purposes, because our re solution is probably insufficient 



to compute the corotation torque (jPaardekooper et al 



20101) • While the interaction with a 



Gl-active disk does not seem to result in a drastically different overall migration rate, the 
direction and magnitude of migration in both simulations is quite variable and seems to 
depend strongly on interactions of the planet with discrete modes in the disk. In t = 10 
planet run, the overall radial migration rate is similar but is even more erratic in direction 
than in the t = planet run. The important departures from laminar disk theory are 
threefold: the more erratic nature of migration, the importance of interactions with discrete 
low-order modes, and stalling of migration near the inner edge of the strongest GI activity. 
The migration is neither a random walk nor monotonic but has a chaotic character due 
to nonlinear interactions. For all these reasons, it will be important in future studies to 
explore the dependence of the planet's behavior on its mass and on its placement relative 
to the site of Gl-eruptions (for t = type runs) and relative to the CRs of low-order GI 
modes (for t = 10 type runs). 

In both simulations, the inward migration appears to halt near 17 AU. This may 
simply be a coincidence, because the overall average migration rates are similar, the starting 
positions are essentially the same, and the calculations have a similar duration. We plan 
to extend simulations like these to much longer times in the future. However, the final 
radial positions of both planets happen to be the location of the inner Lindblad resonance 
(ILR) for a strong m = 2 mode present in both runs. This ILR is associated with a surface 
density enhancement in the disk. Near this enhancement, the surface density gradient may 



be falling steeply enough to 



lalt m igration, a mechanism that can work in a laminar disk 



(jPaardekooper fc Papaloizou 



20091 ). Alternatively, at 17 AU, the level of GI activity is 
quite different interior to the planet {Q > 2) and exterior to the planet {Q ^ 2), and this 
may affect the balance of the Lindblad torques exerted by the disk on the planet. Whether 
the migration really does stall and what the mechanism may be are subjects for further 
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research. 

Summary. The simulations in this Letter show that a planet of only IMj placed near 
the Q-minimum in a O.14M0 disk can hasten the onset of gravitational instabilities and 
change the nature of the dominant modes during the initial burst. Although the magnitude 
of r adial migration (a fe wxlO~^ AU yr^^) is not too different from that computed using 



the 



Tanaka et al. 



(I2OO2I ) isothermal disk formula based on the Lindblad torques, it can 
fluctuate in sign on time scales on the order of the orbit period and can be strongly affected 
by gravitational interactions of the planet with discrete Gl-driven spiral modes, leading to 
phases of rapid radial motion. Planet orbit eccentricities fluctuate between about 0.02 and 
0.10. The simulations also suggest that planet migration may stall for long periods near the 
ILRs of dominant global Gl-modes. 
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Fig. 1. — (Color figure available in the online version) Global amplitudes of nonaxisymmetric 



'or individual cos fmc/)) perturbations. The 



20061 ). Here the integrals are 



density perturbations function of time 

formula for computing is equation (15) of (iBoley et al. 
done only over the disk outside 15 AU to suppress contributions from edge effects in the 
inner disk. The top and bottom panels correspond to the fiducial run (no planet) and the 
t = planet run, respectively. 
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Fig. 2. — Plot of the planet's radial position as a function of time Wp{t) for the two planet 
simulations. The dashed-diamond and solid-star lines (colored red and black in the online 
version) indicate the curves {tu{t)), points (a), and horizontal scales (t) that correspond to 
the t = and t = 10 planet runs, respectively. The symbols are the approximate semi-major 
axes a computed for each to 20 change in azimuth of the planet as explained in the text. 
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Fig. 3. — Plot of the eccentricities e, computed for each to 2(f) change in azimuth of the 
planet as explained in the text, as a function of time. The diamonds and stars (colored 
red and black in the online version) indicate the points (e) and horizontal scales [t] that 
correspond to the t = and t = 10 planet runs, respectively. 



